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Abstract 

We present a new density-functional method of the self-consistent electronic- 
structure calculation which does not exploit any local density approximations 
(LDA). We use the exchange-correlation energy which consists of the exact 
exchange and the correlation energies in the random-phase approximation. 
The functional derivative of the correlation energy with respect to the density 
is obtained within a static approximation. For transition metals, it is shown 
that the correlation potential gives rise to a large contribution which has the 
opposite sign to the exchange potential. Resulting eigenvalue dispersions and 

the magnetic moments are very close to those of LDA's and the experiments. 
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The optimized-effective-potential (OEP) method for the electronic-structure calculations 
was first applied to atoms by Talman and Shadwick who recognized it as a kind of 
restricted minimum search for the Hartree-Fock total energy. In the method, the variational 
space was restricted to the space of the local one-particle potential, which generates all the 
eigenfunctions, or to the space of the corresponding density. From the view of the density- 
functional (DF) theory, their calculations are recognized as the Kohn-Sham (KS) exact 
exchange (EXX) -only DF calculations. In the last years, we have extended their method 
to be applicable for solids [§- 0], where we add the correlation energy in the local density 
approximation (LDA). As was shown there, the results admittedly were not satisfactory 
from the view point of comparison with experiments. Especially for transition metals such 
as iron [Q, the method gave the occupied d bands which were too deep relative to the s bands 
and also gave rise to too large magnetizations. This indicates that the LDA correlation is 
poor when we combine it with EXX, and that the true correlation should give rather large 
contributions cancelling large EXX. 

There are two kinds of way to take into account the correlation energy more precisely. 
One is the multiple-configuration extension which is a widely-used idea in the quantum 
chemistry. The other is based on the DF theory, to utilize the implicit DF for the correlation 
energy based on the KS orbitals, in addition to the EXX energy. Along the second, Garbo 
and Gross tried a new DF for correlation, which applied to atoms and molecules [^. 

In this Letter, we present a new method for solid along the second line. We use an 
exchange-correlation (XC) energy functional, which consists of the EXX energy and a corre- 
lation energy in the random-phase approximation (RPA). The correlation potential, which 
is the functional derivative of the correlation energy with respect to the density, is evaluated 
within a static approximation; the static screened- Coulomb interaction, which is used for 
the evaluation of the RPA correlation, is calculated by use of the product-basis method 
developed by Aryasetiawan and Gunnarson ||^. We denote our method as EXX+RPA in 
the following. 

Formally, the XC energy -Excf?^] as a DF is given by use of the coupling-integral method: 
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where ip denotes the electron field operator, and (...)«, a indicates the vacuum expectation 
value with respect to the system with the given density ^^-(r) and the coupling Ae^. Note 
that the given density is kept unchanged for any A by applying the A-dependent external 
potential. The EXX energy Ey\n] as a functional of the density is defined as the term of 
E^c[n\- As a correlation energy functional i^ci^^], "we only take the RPA contribution included 
in the expression of Eq. (|1|) : 

^^c^'^'^M = ^ /' ^Tr[(l - v^D'Y'v^D' - vy^D'] = ^Tr[log(l - vD') + vD% (2) 
z JO A 2 

occ. unocc. f ^ -i 

D°(ri,r2,c.)=^ Y: C(ri)V^;(r2)^r(r2)^i(ri) , , , , ' ^ 1 ^ , (3) 

i j [uj - eir + el + to uj + tir ~ el - 10 ) 

where the trace is taken for index 1 = rtcr (in Eq.@, we suppress a factor 1/ dt for sim- 

phcity). v\ is defined as V\{1, 2) = \e^5{ti — t2)5^-^(,.J\Yi — r2|, and v denotes f a=i- Formally 

Eq.@ is the same as that for the RPA correlation energy of the homogeneous electron gas 

0] except that E^^^ should be treated as a functional of density riaij) through the Lindhard 

dielectric function , which is spin-diagonal and constructed from the eigenfunctions ipl{v) 

satisfying 

[:^ + r;«(r)-6^R(r) = 0, (4) 

occ. 

^.(r) = ElC(r)r- (5) 

i 

Under the assumption of one-to-one correspondence between no-(r) and V^^(r) (we define 
V^^ so that the chemical potential is included in it), we can treat ^Jlij), and hence E^^^, 
as a functional of ^^-(r). The approximation of Eq.(|^) was suggested by Gross, Dobson and 
Petersilka in Ref. (we omit the term /xc corresponding to the vertex correction). 
The derivative SE^. with respect to 5D^ can be written as 

5E, = '-i:i[W^5D% (6) 
= V,, - V, (7) 
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where Vsc = (1 — vD'^) denotes the dynamical screened Coulomb interaction in RPA. 
We evaluate 6Ec in a static approximation, i.e., we replace PVp(ri,r2,ti — ^2) with 
iyp^=0(ri,r2) X 5{ti - t2), where we define W^=\ri,r2) = JZodtWp{ri,r2,t). This ap- 
proximation is justifiable if the relaxation time of the dynamical screening Vsc, typically the 
plasma oscillation time-scale, is sufficiently shorter than that of the density fiuctuation iD^. 
We know 2D°(ri, r2, ti = ^2) = na{ri)6{ri — r2) — [rio-(ri, r2)]^. Here we use the non-local 
density no-(ri, r2) = J2T^' V^"(ri)'?A^(r2). Then we obtain 6Ec = SEd +6Ec2, where we define 

SEd = ^J2 [ dv,dv2W^=%T,,r2)5iMr^,r2W), (8) 

SE,2 = \T.I dvW'^=\v, v)5n,{v). (9) 
^ 0- 

5Eci and 5Ec2 correspond to the correlated part of the screened exchange and the Coulomb- 
hole terms, respectively (see p. 40 in Ref. [|T3])- As for we can calculate its func- 
tional derivative 5Eci/ 5Vf^{r) from "°(ri, Y2), n^ijci, Y2) and dn^ijci, r2) / 5Vf^ {r) through 
Eq.(H). Then we obtain 5Eci/ 5nfj{Y) by the same inversion method as was used to obtain 
6E^/6ncr{r) (see Eq.(|TTD). If we evaluate 6Ec from Eq.(|]) without the static approxima- 
tion, we obtain the result essentially equivalent to Eq.(27) in Ref. |Tl|, which was used by 



Godby, Shliiter and Sham to discuss the eigenvalues of the density functional theory beyond 
IDA. 

In order to calculate all the related quantities, we exploit the atomic sphere approxima- 
tion (ASA). Any points in the space are denoted by (r, R), where R is the index for atomic 
sphere (AS) and r = (r, 6, 0) is a vector denoting the position in each AS. We consider 
Exc[n^] as a functional of rf^{r,R), where n^^{r,R) denotes the spherically-averaged density. 
Wp^^ are calculated by the product-basis method as 

W^='{v,,R,,V2,R2) = Y.Wr'ihRi,J,R2)HT,)B,ir2), (10) 

where i?j(r)'s form the product basis (see Ref. for notation). In the method of the 
LMTO-ASA |]r2|, the non-local density n^{ri, r2, -R2) is given as a functional of V^^{r, R) 



through the potential parameters and radial basis functions. Therefore 6Eci/6V^ {r, R) can 



be evaluated by use of the derivative chain rule in the same manner used for 0. Now, 
we can calculate the contribution to the correlation potential VJr^(r, R) = 6Eci[n^]/Sn^^{r, R) 
by the inversion equation, 

SVf{r,R) y7o SVf{r,R) 5n%{r',R'y ^ ' 

Eq.(pTD is essentially the same as the one for E^ in Ref. 0. Adding spin-independent 
V'^'^{r,R) = (the spherically-averaged W^^^^r, R,r, R)), we obtain the RPA correlation po- 
tential V^. Since we treat only metals, we can determine constant parts of and V^^ 
uniquely (note that V^^ includes the chemical potential). We have developed a code to 
perform the non-relativistic self-consistent calculation with together with V^. This is 
done by combining two codes, one calculating Wp and the other LMTO-ASA EXX code; 



the former is a part of the GW program |]13l provided by Aryasetiawan, and the latter 
is developed starting from LMTO-4 [p!5| . 

We show results for Cu(fcc), Ni(fcc), Fe(bcc), and Co(fcc) with the lattice constants 
determined by LDA , 6.76, 6.55, 5. 27, and 6.54 a.u. [|l^, respectively. For the evaluation of 
6E^ and SE^, we take all the pairs (i?i,i?2) within the second nearest neighbors and treat 
4s, 4p, and 3d as valence orbitals. For the calculation of W^^^, we use 35-k points in the 
irreducible Brillouin zone, with the Lorenzian broadening 6 = 0.02 Ry. in Eq.(^. We used 
90 optimum product basis -Bj's, where all the core eigenfunctions are used to generate them. 

In Fig.|l], we show the self-consistent V^'^ = V"^ + V^^ + V'^ for Cu, Ni and Fe. Constant 
parts of and are properly fixed for presentation. We find that the large dips of 
and the large differences between up- and down-spin components of are strongly reduced 
by adding . As a result the final V^'^ obtained by EXX+RPA gets much closer to those 



by LDA than it was for obtained by EXX. which corresponds to the screening 



length [|18[, has considerable structures but its contribution to is largely cancelled by 
V^^ in the vicinity of the core region. This can be explained by the fact that is 
sufficiently short-ranged for given bnij) near the core region, which allows us to evaluate 
bEc approximately by use of Vr^"°(ri,ri) instead of Vr^"°(ri,r2) in Eqs.(|-||): By this 
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replacement, we have SEd + 5Ec2 = because of / driD^ijci, r2, ti = ^2) = 0. On the other 
hand, we see a rather large cancellation between and V^^ out of the core region (r ^ 1.0 
a.u.). This is because the main contribution to 5Eci in Eq.(|[) comes from the integral out 
of the range of v^^^, where we expect the behavior of is VF^^''(ri, r2) ~ — l/|ri — r2|. 

The eigenvalue dispersions are shown in Fig.|^ together with those calculated by usual 
LDA, and EXX. Results by EXX+RPA are very close to those obtained by LDA, and are 
much different from the results by EXX. The magnetic moments shown in Table | are also 
in good agreement with the LDA's and experiments. To obtain numerically accurate values 
of the magnetic moments, it is necessary to extrapolate the series of fully converged (about 
the number of k points) results for different 5 towards 5 = 0. This, however, seems too 
time-consuming in the present stage of our unoptimized computer code, and we only show 
the results for 5 = 0.04 Ry. with 35-A; points in Table |I| for reference. 

Let us discuss some possible sources of the difference between our results and those of the 
true DP. These are (1) the product-basis expansion, (2) RPA, and (3) static approximation 
for RPA. As for (1), our test for a model interaction Vsc = exp{—K,r)/r (in this case we 
can obtain accurate results without using the product-basis expansion) indicates that the 
potential Vc is poorer when getting closer to the AS's center. The oscillation of K near the 
nucleus seems to be an artifact due to the method of expansion. However, this fortunately 
little affects on the energy bands and the magnetic moments (we will show the details in the 
subsequent paper). To take into account the effects beyond RPA, (2), we have to treat /xc 
in Ref. 0, which corresponds to the vertex correction. Some parts of the contribution due 
to /xc might be taken into account by LDA-like approximations; one of the simple ways is 



to add the difference between the LDA XC calculated by RPA |17] and the one calculated 



by a more accurate scheme |T9|. We tried the above method but the correction turned out 
to be rather small; the magnetization of Fe enhanced by 0.05 /^b (this value is similar to the 
corresponding LDA case). However, we are not so confident whether the correction is really 
meaningful or not. As for (3), the dynamical effects will be simulated in effect by making 
closer to v (no relaxation limit). It should reduce the magnitude of V^. Therefore 



the position of d band relative to s band should be somehow pushed down for Cu case (See 
Fig.|l]). Contrary, this estimation concerning (3) is opposite in the case of LDA. In the case of 



LDA, we can rather easily show that d bands calculated by LDA in the RPA ||T^ are pushed 
up from those obtained by its static approximation. In conclusion, it seems rather difficult 
to evaluate (3) based on the LDA idea; we have not succeeded in giving any reasonable 
evaluations for the magnitude of errors due to (2) and (3). 

It sometimes have been claimed that the itinerant-magnetism is a sensitive problem 
where the higher-order electron-correlation |^ beyond RPA screening might be essential 



in determining the magnetic moments or the ground states. Our results, however, indicate 
that the essential part determining the magnetic moments might be well accounted for by 
the correlation within the RPA screening, even though the higher-order correlations are 
certainly important for excitations as was discussed in Ref. [1^. Incidentally, as for the the 



higher-order correlations, we as well can extend the OEP method to include other classes of 
correlation diagrams included in Eq.(|I|) and to evaluate their contributions in an ab-initio 
way. For example, the exchange diagrams of order e'^ that is important to cancel the self- 
interaction of order can be included. Further, the ladder types of diagrams that include 
the free energy due to the spin wave should be taken into account in the case of finite 
temperatures. 

Our method determines not only V^^, but also W^^° in a self-consistent manner. W^^° 
calculated from the LDA eigenfunctions (and eigenfunctions) seems reasonable for transition 
metals: For example, the magnetic moments 2.07 /is of Fe, which is calculated with a fixed 
Wp^^ calculated from the LDA eigenvalues, is very close to the self-consistent value 2.04 
/is in Table III However, such results obtained by the LDA eigenvalues are not expected to 



be sufficient for solids like transition- metal oxides. In Ref. |T^ for NiO, it is shown that the 
GW calculation with Wp constructed from the LDA eigenvalues gives poor results because 
of the too small a LDA band gap, and that the self-consistency of GW is essentially needed. 
The fully self-consist calculation of GW seems very hard mainly due to the programing 
difficulties and the lack in computational power. From such a point of views, our self- 



consistent calculation might be a good substitution of the self-consistent GW. 

We have not developed the code to calculate Ef^^ itself yet. We only treated its func- 
tional derivative dE^^^/Sn^, and the static approximation is not applicable for the total 
energy calculation. To our knowledge, none have done the ab-initio calculation of Ef'^^. It 
requires to treat unoccupied states precisely. It might be a httle difficult in LMTO-ASA 
scheme. However, we have to try such calculations for determining lattice constants and 
other cohesive properties. 

I thank F. Aryasetiawan for his GW program and for stimulating discussions, and M. van 
Schilfgaarde, T. A. Paxton, O. Jepsen, and O. K. Andersen for their TB-LMTO program. I 
also thank H. Akai for discussions and for critical reading of the manuscript. 
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FIGURES 

FIG. 1. Exact exchange potential and RPA correlation potential = V^^ + V^"^ for Cu, 
Ni and Fe. The LDA XC potential |17| as references are calculated for the density determined by 



the self-consistent calculation of EXX+RPA. 

FIG. 2. Energy dispersion curves calculated by three different types of self-consistent calcu- 
lations, LDA, EXX, and EXX-I-RPA. EXX denotes the calculation which is performed with LDA 
correlation [^] potential in addition to the EXX potential. 
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TABLES 

TABLE L Calculated spin magnetic moments (/xb)- For EXX+RPA, we show the values 
obtained with 6 = 0.02 Ry, and 6 = 0.04 Ry (parenthesis), see Eq.(|3|). The lattice constants used 
in the calculations are somehow smaller than the experimental ones (see text). Experimental spin 
magnetic moments are deduced from the experimental values of the saturated magnetization and 
g- values 

LDA EXX EXX+RPA Expt. 

Fe 2.13 3.27 2.04 (2.07) 2.12 

Co 1.54 2.29 1.51 (1.52) 1.59 

Ni 0.58 0.68 0.57 (0.58) 0.56 
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